Mon. Not. R. Astron. Soc. 000, [T]{T3] (0000) Printed 4 March 2013 (MN M£X style file v2.2) 



Crab Nebula gamma-ray flares as relativistic reconnection 
mini jets 

E. Clausen-Brown 1 * and M. Lyutikov 1 ! 

1 Department of Physics, Purdue University, 525 Northwestern Avenue, West Lafayette, IN 47907-2036, USA 



Accepted/Received 



ABSTRACT 

The unusually short durations, high luminosities, and high photon energies of the 
Crab Nebula gamma-ray flares require relativistic bulk motion of the emitting plasma. 
We explain the Crab flares as the result of randomly oriented relativistic "minijets" 
originating from reconnection events in a magnetically dominated plasma. We develop 
a statistical model of the emission from Doppler boosted reconnection minijets and find 
analytical expressions for the moments of the resulting nebula light curve (e.g. time 
average, variance, skewness). The light curve has a flat power spectrum that transitions 
at short timescales to a decreasing power-law of index 2. The flux distribution from 
minijets follows a decreasing power-law of index ~ 1, implying the average flux from 
flares is dominated by bright rare events. The predictions for the flares' statistics can 
be tested against forthcoming observations. We find the observed flare spectral energy 
distributions (SEDs) have several notable features: A hard power-law index of p < 1 
for accelerated particles that is expected in various reconnection models, including 
some evidence of a pile-up near the radiation reaction limit. Also, the photon energy 
at which the SED peaks is higher than that implied by the synchrotron radiation 
reaction limit, indicating the flare emission regions' Doppler factors are ;> few. We 
conclude that magnetic reconnection can be an important, if not dominant, mechanism 
of particle acceleration within the nebula. 

Key words: magnetohydrodynamics (MHD) - magnetic reconnection - radiation 
mechanisms: non-thermal - pulsars: general - ISM: individual: Crab Nebula - ISM: 
jets and outflows - supernova remnants 



1 INTRODUCTION 

The presumed constancy of the high energy Crab nebula 
emission has surprisingly been shown to be false by multiple 
day- to week-long flares, presenting a challenge to standard 
pulsar wind models (; Kennel fc Coroniti|1984 1. During these 
events, the Crab Nebula gamma-ray flux above 100 MeV 
exceeded its average value by a factor of several or higher 



(Abdo et al. 2011 


Tavani et al. 2011 Buehler et al. 2012 


|Striani et al.||2011 


, while in other energy bands nothing 



unusual was observed (e.g. Abdo et al. 2011 Tavani et al 



2011 and references therein). Additionally, sub-flare vari- 
ability timescales of ~ 10 hours has been observed ( |Balbo| 
et al |2011||Buehler et al.|2012 i. 

There are two interesting observational facts related to 
the gamma-ray flares. First, is their unusually short duration 
of a few days. This time scale, on the one hand, is millions 
of times longer than the period of the pulsar, yet on the 
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other hand it is hundreds of times shorter than the nebula's 
dynamical timescale of ~ few years. We consider it unlikely 
that the flare is related to the changing plasma properties 
within the pulsar's light cylinder, both due to the extremely 
large separation of temporal scales and due to the fact that 
no changes in the radio pulsar timing properties were seen 
during the flare (Espinoza et al. 2010[ ). Thus, we associate 



the duration of the flare with stochastically changing prop- 
erties of plasma within the nebula. Second, the flaring be- 
havior consists of apparently isolated, intermittent events 
that are dominated by bright rare flares. 

Two contrasting models of these flares can be envi- 
sioned. First, a flare can be due to large scale changes in 
the steady nebular flow, amplified by the effects of relativis- 
tic beaming. Operating within this model, Komissa rov fc| 



Lyutikov (2011 \ and Lyutikov et al. (2011 1 place the flaring 



location in the downstream region of an oblique shock. The 
post-shock flow of an oblique shock can be highly relativis- 
tic, and thus Doppler boosted, with a bulk Lorentz factor 
of r ~ <f>~ 1 , where <f> s is the angle between the upstream 
velocity and the shock plane. The short timescale can be 
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attributed to the shock normal changing direction (perhaps 
due to sausage or kink instabilities, or a corrugation per- 
turbation), causing the post-shock velocity to sweep across 
the line of sight and create a short flare in Doppler boosted 
emission (the lighthouse effect). 

Second, the flare can be due to a highly localized emis- 
sion region, or blob, so that the flare observables determine 
the intrinsic properties of the emission region. The natural 
flaring mechanism in this category is relativistic magnetic 
reconnection — the focus of this work — which has been in- 
voked by Crab Nebula flare models ( Bednarek fc Idec|2011 



Cerutti et al. 2012 Uzdensky et al. 2011 1 and fast flaring 



models in gamma-ray bursts (GRBs) and active galactic nu- 



clei (AGN, 


Lyutikov & Blackman|2001| |Lyutikov & Bland- 


ford|2003| |1 


^yutikovpOOei |Giannios et al.|2009[|2010| |McK- 


inney & Uzdensky 2010). 



To see why reconnection is a natural flaring mechanism 
in PWNe, consider that the energy budget of the pulsar wind 
is set by the spin-down power, P sp i n , of the pulsar. In the 
standard model, P sp m is a smooth, monotonically decreas- 
ing function of time, implying that the nebula emission will 
track this smooth decline. In contrast, the emission from 
relativistic outflows in systems such as GRBs and AGNs 
are probably related to an irregular accretion process. Their 
variable accretion rates and other accretion disk instabili- 
ties may produce unsteady outflows with internal shocks, 
providing a viable alternative to reconnection as a flaring 
mechanism in the outflow (e.g. |Rees fe Meszaros|1994[|Spada| 
et al.|2 001). In pulsar outflows, however, the internal shock 
scenario is more difficult to realize given the regular behav- 
ior of P sp i n . Instead, it is more likely that day-scale flaring 
is due to the build up and eventual release of free energy 
somewhere in the PWrsQ Magnetic reconnection is a natu- 
ral candidate for this process as it is the primary mechanism 
by which magnetic free energy may be suddenly converted 
into particle kinetic energy, which occurs via the topological 
rearrangement of magnetic field lines. Furthermore, in the 
closest space laboratory available, the solar system, this is 
precisely what is observed, where intermittent flaring and 
the production of non-thermal particles is regularly associ- 
ated with magnetic reconnection ( Priest fc Forbes|2000 1. 

In this work we locate the reconnecting plasma at mul- 
tiple sites in the nebula, and presume the reconnection 
plasma's Alfven velocity approaches the speed of light. The 



reconnection outflow speeds are then relativistic (Lyutikov 



& Uzdensky 2003 Lyubarsky 2005 1 and behave as "mini- 



jets," a model that has been used to overcome the gamma- 
ray opacity problem in the context of gamma-ray bursts and 
active galactic nuclei ( Lyutikov & Blandford 2003 Lyutikov 



2006 Giannios et al. 112009 1. In an approach similar to ours, 



Yuan et al. (20111 invoke a series of relativistically mov- 



1 The wisps in the nebula that exhibit variability timescales of 
~ months, are not due to an unsteady wind, but are probably re- 
lated to the ion-affected structure of the wind termination shock 
(|Gallant fe Arons|1994} or a synchrotron instability ( Hester et al.| 
|2002| |, However, the recent discovery of the several year ~ 10% 
decrease in the nebular X-ray emission (Wil son-Hodge et al.|2011| 
is puzzling in light of the Crab's smooth spin-down rate, but the 
unknown mechanism for this emission decrease is probably unre- 
lated to day-long gamma-ray flares, since its timescale is several 
orders of magnitude longer. 



ing "knots" ("minijets" in our terminology), with statistical 
properties that are compared to the Crab Nebula light curve 
via Monte Carlo simulations (they do not identify the knots 
with reconnection). Our model differs from that of Yuan 



et al. (20111 in that it is analytical, and it presumes that 
relativistic beaming controls the observed statistical prop- 
erties of the reconnection minijets. 

This paper is organized as follows. In j|2]we first focus 
on constraints to the emission regions' magnetic field derived 
from adiabatic expansion and synchrotron cooling, and com- 
pare these to Crab Nebula magnetic field estimates based on 
cquipartition and spectral modeling. We include the possi- 
bility that the emission region has a bulk relativistic velocity 
with a Lorentz factor of V and associated Doppler factor of 
5 = (r — \/F 2 — 1 cos^oi,)^ 1 , where 8 t is the angle between 
the blob velocity and the line of sight. Due to the high lumi- 
nosity and short variability of the flare, we assume o b ~ 0, 
in which case 8 ~ 2T. We then speculate in ij3]that the Crab 
Nebula flare originates from a relativistic outflow caused by 
magnetic reconnection in a high-cr plasma (a = Poynting 
flux/particle flux, Kennel & Coroniti 1984 1. ^4] focuses on 
the emission region spectral energy distribution (SED) in the 
context of the synchrotron radiation reaction cut-off and the 
functional form of the observed SEDs . In Sj5j we construct 
an analytical statistical model which describes the statistics 
of both individual flares and the entire light curve's sta- 
tistical moments (e.g. average and variance) and its power 
spectrum. We provide a summary of our conclusions in fj6j 



2 EMISSION REGION PARAMETER 
ESTIMATES 

In this section the magnetic field of the emission region, or 
blob, is constrained under two separate assumptions: that 
(i) adiabatic cooling or (ii) synchrotron cooling controls the 
flare duration. The relevant observed flare parameters are 
the flare duration, r = 10 r 5 . 5 s; typical photon energy, 
e p = lOOeioo MeV; the blob's isotropic equivalent luminos- 
ity, Li SO = 10 36,6 L36.6 ergs s _1 ; and blob magnetic field, 



B — 0.3-B cq G; where the fiducial flare parameter values are 
based on the Fermi/LAT data of the September 2010 flare 
(Abdo et al. 20111 and the nebula equipartition magnetic 
field ( |Trimble|1983| ). 

Assuming adiabatic cooling limits the flare duration, a 
lower limit on the blob frame magnetic field can be made by 
examining the energy content of the blob. We parametrize 
the total magnetic energy of the blob as being some sig- 
nificant fraction /m of the blob's total energy content, 
JmE' = B' 2 R' Z /6. Assuming the blob is causally connected 
and is moving directly along the line of sight, the blob ra- 
dius is R' ^ cFt. Of the total energy content of the blob, 
E = f^TR' 3 B n /&, if a fraction f rad is radiated away dur- 
ing the flare, then the total radiated power is L — fradE/r. 
L is related to the isotropic equivalent luminosity, L 
assuming that all of L is beamed into a solid angle of 7rr 
which leads to the relation, Lt so = 4r 2 L. Thus, the isotropic 
equivalent luminosity can be expressed as 



iso, by 

2 



-lp3 



(1) 



In order to use equation |TJ to constrain the magnetic field, 
R' is estimated, as before, using the light crossing time r > 
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R'(Tc) 1 . This constraint may be interpreted as the blob 
slow down time (see, e.g., Giannios et al.||2009) or as the 



adiabatic cooling time. That is, if the blob begins with an 
initial radius R' and expands at near light speed, R' t = R' + 
ct' then the blob light crossing time, R' /c, is the time the 
blob takes to double in radius and adiabatically cool by a 
significant amount. Solving equation |TJ for B' and using 
the inequality R' err allows us to set a lower limit on B': 



B' > 1.5(/M// rad ) 1/2 L^ 2 6 r^r- 



mG. 



(2) 



If we assume 10% of the blob's energy content is radiated 
away (f ra d = 0.1) and the blob is magnetically dominated 
(/m ~ 1), then B' > 4.7 mG. Typical PWN models suggest 
the pulsar wind begins magnetically dominated, and some- 
how transitions to being cold and particle dominated near 
the termination shock so that the magnetization parameter 
is cr = B' 2 (8irp'c 2 )~ 1 = 3x 10 -3 , where p is the mass density 
( Kennel & Coroniti 1984 1 . However, exactly how a evolves 
from a S> 1 near the pulsar to a < 1 at the termination 
shock is difficult to explain, provoking competing models 
that suggest cr > 1 at the shock ( Begelman||1998 Lyutikov 
2010). Thus, near the inferred location ot the reverse shock 
the blob could be magnetically dominated so that /jy/ ~ 1 
(in a cold plasma, }m = a/(l + a)) and B' 3> the standard 
nebula value (0.3 mG), or in the canonical model where the 
plasma is not magnetically dominated, fu may be small and 
B' within the measured range of ~ 0.1 mG. 

The synchrotron radiative cooling timescale sets a more 
robust lower limit on B' because it avoids the uncertainty 
in f r and /m that appear in equation |2|. The synchrotron 
cooling time, r c = 5T0 8 (2r7_B' 2 ) -1 s, can be found given the 
typical synchrotron photon energy of e p = 2F'y 2 heB' / (m e c), 
where 7 is the lepton Lorentz factor, h is Planck's constant, 
e is the elementary charge, and m e is the lepton mass. These 
two expressions combined give a synchrotron cooling time of 



2i f - 1 /2 B '-3/2 r - 

Zle 100 -"cq 1 



1/2 



days. 



(3) 



If a significant fraction of the blob energy is dissipated via 
synchrotron radiation in the flare, then r ^ r c and a lower 
limit can be set on the blob frame magnetic field: 



B' > 0.97e- 



-1/3 -2/3p-l/3 



mG. 



(4) 



Note that this lower limit (which is similar to that found in 



Bednarek & Idee 20111 exceeds the observed values of the 



nebular magnetic field (0.1 to 0.3 mG), unless relativistic 
motion is invoked. 

The upshot of these parameter estimates is that the 
magnetic field of the emission region appears to be signifi- 
cantly higher than both nebula equipartition field estimates 



(0.3 mG, Trimble 1983 I and measurements based on mod- 



eling of the nebula SED (0.1 mG, Abdo et al. 20101. This 



suggests that the flaring either occurred in a region of high 
magnetic field in the nebula, or that the emission region is 
moving toward Earth with V ^> few, which reduces the above 
magnetic field estimates. 



3 MAGNETIC RECONNECTION? 

Magnetic reconnection provides a natural explanation of 
the implied relativistic motion discussed above, the intrinsic 



short timescales, and the flares' intermittency. Reconnec- 
tion is a process in which the magnetic energy of a localized 
region, a current sheet, is suddenly converted to random 
particle energy, and bulk relativistic motion (for studies on 
relativistic reconnection, see |Blackman fc F ield 1994; |Lyu-| 



tikov & Uzdensky||2003| |Lyubarsky||2005l lUzdenskyl 2011 1 



With regard to PWNe, reconnection has already been stud- 
ied as a possible resolution of the well known cr-problem 



(Coroniti 1990 Lyubarsky & Kirk 20011. In the canonical 



PWN model ( |Rees fc Gunn|1974||Kennel fc Coroniti|1984| ), 
the magnetization parameter a is high only for radii that 
are well within the wind termination shock, r s , and reduces 
to a ~ 10~ 3 to 10 -2 near r s . To explain how a is reduced 
so drastically as the plasma propagates out to r 3 , reconnec- 
tion in a striped wind has been invoked as the mechanism 
by which magnetic energy is transferred to particle energy, 
thereby reducing cr. In a challenge to the canonical PWN 
model, Begelman ( 1998j ) argues that the toroidally domi- 
nated large-scale nebular magnetic field is subject to the 
m — 1 kink mode instability near the inferred location of r s , 
causing the nebular field to have coherence lengths on the 
order of r„ instead of the size of the radio nebula as it is in 
canonical models. This obviates the need for such low cr at 
r s and may cause reconnection throughout the nebula. In 



a similar vein, Lyutikov (20101 proposes a model in which 



reconnection occurs primarily along the rotation axis and 
equatorial region well beyond the light cylinder, thus quali- 
tatively reproducing the jet/equatorial wisp morphology of 
the nebula. 

In our simple model, we propose that multiple recon- 
nection sites are located in a region of the inner nebula with 
a magnetization parameter of order a ~ few, thereby ac- 



celerating mildly relativistic outflows (Lyutikov & Uzden- 
|sk71|2003"l |Lyubarsky|[2005l ). We assume the nebular mag- 
netic field is similar to the reconnection minijet's emission 
region magnetic field, B nc t ~ B'. The observational nebu- 
lar magnetic field estimates range from -B ne b ~ 0.14 to 0.3 
mG (|Trimble||1983| |Abdo et al. |2010[ ). These values of B ncb 



are clearly less than the lower limit set by the cooling equa- 
tion Q of B' > 1 mG unless we impose relativistic motion. 
Solving for the Lorentz factor in equation Q, we find 



T > 34£ 



10015. 5 J 



(5) 



Such a high Lorentz factor may not be required if the mag- 
netic field is higher near the reconnection region. In split 
monopole models of PWNe, the toroidally dominated mag- 
netic field (B<f, oc sin 8, where 9 is the polar angle) is high- 
est in the equatorial regions ( |Michel|1973}|Bogovalov|1999[ ). 
Thus, if reconnection occurs in the equatorial current sheet 
or equatorial striped wind, then the magnetic field will in- 
deed be higher then the nebula average of 0.14 — 0.3 mG. In 
MHD simulations of the Crab pulsar winds, some parts of 



the nebula magnetic field reach 0.45 to 0.6 mG (Komissarov 
& Lyubarsky 2004 ; Komissarov fc Lyutikov|2011 1 which re- 



quires that the Lorentz factor be at least 4 to 10 according 
equation jjjk. We note that it is also possible there is no 
relativistic motion, and the flaring region simply has a high 
magnetic field of ;> 1 mG, a value that is indeed found in 
bright localized regions of the nebula (e.g. the inner "wisps" 
and "knots" of the nebula, Hester et al.|1995 1 . 
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4 MINUET SED 

The spectral energy distribution (SED) produced in a flaring 
region may be significantly influenced by the proximity of 
the emitting particles to the synchrotron radiation reaction 
limit implied by MHD such that the electron distribution is 
either mono-energetic or a power-law with an abrupt cut- 
off. If a mono-energetic electron distribution is not formed 
and magnetic reconnection is what causes the flares, then 
reconnection may leave an imprint on the power-law index 
of the accelerated particles in the form of a hard electron 
distribution. Here we discuss these features in relation to the 
observed SEDs from the April 2011 flare and suggest these 
observations point toward relativistic reconnection minijets. 

The best Crab Nebula flare SED observations to date 
are 11 SEDs taken during the April 2011 flare by the 
Fermi/LAT team ( |Buehler et al.|2012[ ). [Buehler et al.| ( |2012 l 
fit the SEDs with an empirical function of the form 

N € = Ae ,F exp (— e/ e c ), (6) 

where e represents photon energy, and different values for A 
and e c were used for each SED fit. The parameter 7f was as- 
sumed to be constant for all of the SEDs, and its best fit was 
found to be 7f = 1.27±0.12. Significantly, the observed inte- 
grated flux above 100 MeV, F, correlates with the flare SEDs 
peak energy, e pcak , such that d log F/d log e pca k = 3.42±0.85 



as expected in simple models of Doppler beaming (Lind & 



Blandford 19851. In the context of the above observations, 



we discuss the development of an effectively mono-energetic 
electron distribution near the radiation reaction limit from a 
hard electron distribution of the type expected in magnetic 
reconnection. 

Mono-energetic distributions have already been exam- 
ined in the context of the Crab flares in |Uzdensky et al.] 
( |2011[ ) and |Cerutti et al.| ( |2012[ ) , in which a specific model 
of reconnection causes electrostatic particle acceleration that 
produces an approximate mono-energetic electron distribu- 
tion near the energy associated with the electrostatic poten- 
tial drop. Unlike the model we discuss here, the |Uzdensky] 
et al. (20111 model avoids the MHD limit on particle en- 



ergy and displays no bulk relativistic motion by invoking a 
specific geometry of the reconnection region. In this section, 
we do not make such geometrical requirements, instead we 
focus on a general discussion of how the radiation reaction 
limit can affect flare SEDs. 

A critical feature in the SED of the Crab Nebula flares 
that has already received attention is the synchrotron emis- 
sion that is above the classical synchrotron limit of ~ 10 2 
MeV ( |de Jager et al.|1996|[Lyutikov|2010"l|Abdo et al.|2011 1. 
Importantly, the electron distribution function may display 
a power-law with an excess, or pile-up, of particles near the 
synchrotron limit, in which case the emitting particles will 
display a SED that is close to the single-particle synchrotron 
SED. This limit comes from ideal MHD, in which the electric 
field, E, is less than the magnetic field such that E = r)B, 
where < r\ < 1. To illustrate these points, we briefly ex- 
amine the particle acceleration process. First, we assume a 
particle's energy is approximately described by 
dj eE 
dt m e t 

= rjuiB — PsJ 2 
[/(mlc 5 ), E = TjB 



where f3 s = 2/3e 



(7) 

and ujb ~ eB/m e c. 
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Figure 1. Top plot. A solution for N(-y,t) (see eqn. [9] in text) 
is shown for t = 10t acc and 16t acc , with different normalizations 
to ease visual comparison. The SEDs corresponding to N (7, t) at 
t = 10t a cc and I6t acc are in the bottom plot. For t = 10t aC c, the 
solution is a pure power-law of N oc *y~ 11 ; for t = I6t acc , a pile- 
up develops short of 7 = 7 rat j. Bottom plot. Plotted here are four 
different flare SEDs added onto the average nebula SED : equation 
(JgJ with the parameters reported in |Buehler et al.](]2012| , an SED 
derived from a power-law electron distribution of N = K e y~ 1 ' 1 d'y 
and N = for 7 > 7 max (i.e. the t = 10i aC c solution to eqn[9|, 
the SED derived from the pile-up solution to equation |9| for 
t = lGt acc , and a single-particle synchrotron SED. For visual 
comparison, all three SEDs are adjusted so they have the same 
maximum and peak energy, otherwise the t = 10t acc solution 
would peak at a lower energy than the solutions for t = ldt acc 
and t — > 00 (the mono-energetic solution). 



The first term describes particle acceleration by the electric 
field (where E is only an approximation of the true electric 
field since the accelerating particle velocity is not always 
parallel to the electric field) and the second term describes 
synchrotron losses. We expect the initial particle popula- 
tion to be accelerated to higher energies until drf/dt = 0, 
where the synchrotron energy losses equal the acceleration 
rate such that rjWB = PsJ 2 - If ideal MHD holds in the ac- 
celeration region, then rj < 1, which in turn implies the 
existence of a maximum possible Lorentz factor allowed by 
the synchrotron radiation backreaction for a given magnetic 
field value: 



7rad = 



2e 3 Bi 



= 5 x 10V /2 i?] 



1/2 

,cq • 



(8) 



To quantify what we mean by an electron distribution 
function that "piles-up" near 7 ra( j (e.g. Schlickciscr 19841, we 



examine time evolution of the electron distribution function 
N(t,'y)d'y, which describes the number of electrons with a 
Lorentz factor between 7 and 7 + d7 at time t. The time 
evolution of AT (7, t) in our toy model is described by the 
equation (Kirk et al.||1998 1 



ON d 
dt c*7 



7 



- M N 



+ 7 A ^ = Q%-7o), (9) 
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Figure 2. The Fermi/LAT data from the most energetic part of 
the April 2011 flare are shown (figure 6, panel 7 of |Buehler et all] 
|2012[ | with the corresponding best fit curve from equation ||6| and 
SEDs from three different electron energy distributions: a p = 1 
power-law with an abrupt cut-off, the same for p = 2, and a mono- 
energetic electron distribution. All of the curves are summed with 
the constant nebula SED component used in Buchlcr c t a.1. | ( J2012| | . 
The shaded area represents the one-rx error region (see text for 
more details). 



where the acceleration time t acc and the escape time t esc are 
constants, N(y, 0) = 0, Q is the injection rate of particles 
with Lorentz factors of 70, and S is the Dirac-delta. Equa- 
tion |9| describes the time evolution of non-thermal parti- 
cles subject to mono-energetic injection, particle accelera- 
tion and escape, and synchrotron radiation. The solution to 
equation (19 



Kirk et al. 



subject to the above conditions can be found in 
1998 1. In this solution, the electron distribution, 
N, becomes mono-energetic as t — ¥ 00 (i.e. N ~ 5(y ~ 7rad)) 
if the particle injection process ceases and the acceleration 
process continues, since all the injected electrons will simply 
be accelerated up to 7 ra d- Even if the injection process con- 
tinues during the acceleration process, then the distribution 
can still develop a "pile-up" just below 7 ra d, effectively be- 
coming mono-energetic as discussed below (see also fig. [I]). 
However, such a pile-up distribution only occurs depending 
on the details of the injected spectrum of particles and par- 
ticle escape. In the particular case described above, a pile-up 
only occurs if t esc > t acc (see also Schlickei ser|1984[ ). 

We have plotted the solution to equation Q at two 
time-slices, t — lQt acc and t = 16t acc , in figure [I] for the pa- 
rameter choices t esc = 10t aC c and 7 rac i — 4x 10 .As figure [l] 
illustrates, a pile-up does eventually accumulate just short 
of 7 = 7 r ad, eventually evolving to a mono-energetic dis- 
tribution that produces a single-particle synchrotron SED 
as shown in the bottom panel. Visual inspection of figure [I] 
(bottom panel) makes it clear that the power-law distribu- 
tion of N oc 7 _1 compares more favorably than the mono- 
energetic distribution with the best fit of the April 2011 flare 



SEDs repo rted in|Buehler et al.| i ]2012 |. 

Unlike Buehler et al. (2012 I, we examine only the most 



luminous observed SED from the April 2011 flare (figure 6, 
panel 7 in Buehler et al.|2012 |, since over the course of the 
nine-day flare, the functional form of the SED may have 
significantly changed, and this SED is least affected by the 
background nebula emission. As shown in figure [2] we fit the 



SED to equation (Jm, which yields jf — 1.08 ± 0.16 with a 
reduced \ 2 of 0.35. The shaded area in figure |2j) is the "one- 
cr error region." It represents a subset of curves generated 
by equation jSJ, each of whose parameters lie within a one-er 
interva|^]of the best fit parameters. 

Equation (JsJ) is a useful SED fitting function for the rel- 
evant energy range because it has the same functional form 
(to within a few percent) as SEDs derived from three dif- 
ferent electron distributions: power-laws distributions with 
abrupt cut-offs, power-law distributions with pile-ups (as 
in fig. [T|, and mono-energetic distributions. Note, however, 
that the functional form of SEDs derived from power-law 
distributions that are near the tF t peak depend sensitively 
on the form of the electron distribution cut-off. For expo- 
nential cut-offs to the electron distribution function, the de- 
rived SED is broader than SEDs derived from power-laws 
with abrupt cut-offs that are approximated by equation Q. 
Here we presume the radiation reaction limit is the mecha- 
nism whereby either an abrupt cut-off or a pile-up develops 
in the electron energy distribution. Also, near the eF e peak, 
the 7f parameter cannot be interpreted in the usual way as 
being equal to (p + l)/2 for power-law distributions, or 2/3 
for mono-energetic distributions. For example, SEDs corre- 
sponding to p = 1 and p = 2 are well approximated by 
7f ~ 1.30 and 1.56 respectively, while a mono-energetic dis- 
tribution has jf ss 0.7 (e.g. Mchosc 1980); for the purpose 



of comparison, the best fits of these distributions (which 
only vary normalization and cut-off energy since p is fixed) 
are plotted in figure [2] next to the best fit of equation iEJ. 

The best fit value of 7f = 1.08T0.16 lies in between that 
expected from a mono-energetic distribution and the p — 1 
distribution, while distributions for which p ^ 2 are unlikely 
to have produced the observed SED. The precise p corre- 
sponding to the best fit of -y F = 1.08 ± 0.16 is p ~ -0.2±°- 9 , 
a large range of values due to the calculated SED's weak 
dependence on p for p < for the energy range in ques- 
tion. The best fit values for 7f may imply the development 
of a pile-up near 7 ra d from a hard power-law electron dis- 
tribution of p ~ 1 as seen in figure [l] and discussed in, for 
example, Schlickeiser ( 1984 1. A pile- up pushes ap = 1 SED's 



7f value down from jf ~ 1-30 to that of the mono-energetic 
distribution, "/p = 0.7, as illustrated in figure [l] We verify 
that pile-up scenario could explain the data through the toy 
model we develop from equation ([9]) by setting 70 = 10 4 , 
7r a d = 4 x 10 9 , and t CBC 3> t aC c, after which we find the 
best fit 7f can be reproduced by the resulting pile-up that 
develops when t is between ~ 13i acc and ~ 18ta.ee ■ Such a 
non-steady state pile-up scenario has implications for the 
time evolution of the flare SED that could in principle be 



applied to the other SEDs reported in Buehler et al. (2012 1. 
However, we do not pursue these implications here due to 
the other SED's higher contamination by the background 
nebula. Thus, a wide variety of electron distributions could 
have produced the observed SED shown in figure [2] but we 
rule out standard shock acceleration scenarios with p > 2 



2 The error analysis was done via Monte Carlo simulation. We 
assume Gaussian errors in the measured eF e values, then ran- 
domly generate 10 4 synthetic SEDs, and fit each synthetic SED 
to equation (J6J. The resulting spread about the best fit value for 
each parameter in equation (JsJ allowed the calculation of each 
parameter's standard deviation, or a. 
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Figure 3. Shown are Doppler boosted single-particle synchrotron 
SEDs (dotted lines) for 5 = 1 and 5 = 3. The normalization of 
the intrinsic SED and the critical frequency values are fixed so 
that the 5 = 3 SED displays the maximum reported eF e value 
and correct peak energy during the April 2011 fla re: (e-F e ) max ~ 
4 X 10~ 9 ergs cm -2 s* 1 and e pcak = 375 MeV (Buehler et al. 
2012) ■ Th e average nebula SED (thick line) reported in |Buehler et 
al.| | |2012] | is summed with the flare SEDs to produce a combined 
SED shown by the dashed line. Note that the difference between 
5=1 and 5 = 3 SEDs is enough to render the former nearly 
unobservable compared to the average nebula emission. 



and suggest it was most likely a hard distribution of p <J 1, 
possibly with a pile-up near the radiation reaction limit. 

If the electron distribution is a hard power-law of iV oc 
7 , where N = for 7 > 7 ra d, the observer frame peak 
energy in the eF € representation is 



^pcak 



0.6?ia; c (7 r ad) ~ 1.357/ 



hm e c 3 



IOO77 MeV, (10) 



where the critical synchrotron frequency is ui c = 
3/2y 2 eB±/m e c, and the location of the peak energy, 
0.67iw c (7 ra d), was determined numerically. The observed 
peak energy of a Doppler boosted SED is e p0 ak = <5e peak , 
implying that any observed peak energy above 100 MeV 
must have a minimum Doppler factor of 5 min = e P cak/e p eak- 
Thus, for the April 2011 flare, 5 min > 375/100 ~ 4. 

However, if the theoretical criteria discussed above are 
met for the electron distribution to become mono-energetic, 
then the SED will be effectively described by the single- 
particle synchrotron SED, whose form in the low energy 
limit is eF e oc e 4//3 , and in the high energy limit, eF e oc 
e 3,/2 exp (— e/e c ). Since the eF e mono-energetic synchrotron 
spectrum peaks at the angular frequency of cj pC ak ~ 1.3u c , 
this implies the tF t peak from emitting particles with 
Lorentz factors of 7 ra d is located at a photon energy of 

YlTfl r 3 

4cak « 1.3fiw e (7«d) « 2.97/ c — « 2007/ MeV. (11) 

Thus, the April 2011 flare implies a minimum Doppler fac- 
tor of 5 min > 2. Note that this is a lower estimate of 5 min 
compared to the estimate made assuming a hard power-law 
electron distribution. 

To illustrate the importance of Doppler beaming, we 
have plotted in figure |3] the same intrinsic single-particle 
synchrotron SED with two different Doppler factors, 5 = 1 
and 5 = 3. The different Doppler factors affect the in- 
trinsic SED's photon energies, e = Se', and normalization, 



eF e (e) = 5 A e'F' t ,(e') ( |Lind fc Blandford|l985| ). We adjust the 
normalization of the intrinsic SED to ensure that the 5 = 3 
flare has a maximum at eF e ~4x 10 -9 ergs cm -2 s _1 , as 
observed. The same intrinsic normalization is used in the 
5 = 1 flare. The constant Crab Nebula SED as described in 



Buehler et al. (20121 is plotted as well. Note that by using 
Doppler factors consistent with the SED cut-off observed 
during the most luminous part of the April 2011 flare (S ^ 
few), the unbeamed flare does not significantly increase the 
high energy flux of the average nebula. 

The non-detection of significant flaring by X-ray tele- 
scopes places a further constraint on the flare SEDs. If 
the flare SED comes from a mono-energetic electron dis- 
tribution, then in the Chandra and XMM-Newton energy 
bands, which we define here as extending from e m i n = 0.1 
keV to e max = 10 keV, the SED goes as eF e oc e 4/3 . For 
the single-particle synchrotron SED discussed above, where 
(eFe)max ~4x 10~ 9 ergs cm" 2 s" 1 and e pca k = 375 MeV: 

,T0 keV 

F x = A e 1/3 de w 3.24 x 10~ 14 ergs cm" 2 s" 1 

J 0.1 keV 
,10 keV 

N x = B I e~ 2/3 de w 6.39 x 10" 6 cm" 2 s" 1 , (12) 

JOT keV 

where Fx is the integrated X-ray flux and Nx is the total 
photon flux. These values are much smaller than the spa- 
tially integrated Crab Nebula X-ray energy flux of ~ 10~ 7 
ergs cm -2 s _1 and photon flux of ~ 100 cm -2 s _1 (Kirsch 



et al.||2005 1. Even the photon flux of the Crab nebula in a 
Chandra resolution element of ~ 1 arcsec 2 , which is ~ 10~ 2 
cm -2 s , is well above the photon flux we predict in the X- 
ray band. If the flare emission were from a p = 1 power-law 
electron distribution, then the estimates in equation ( 12 1 
increase by a factor of (7rad/7min) 2 '' 3 , assuming the lowest 
energy range of the SED goes as F e oc e 1//3 . Such a power- 
law could extend down to 7 m i n ~ 10~ 5 7 ra d before the flare 
were comparable to the Crab Nebula flux in one Chandra 
resolution element. Hence, with a 7 m i n as low as ~ 10 4 , it 
is possible to explain the non-detection of the flaring events 
by X-ray telescopes. 

Our above discussion of the flaring SEDs has two im- 
plications: (i) for hard electron distributions, the MHD ra- 
diation reaction limit can lead to the formation of a pile-up 
electron distribution that is effectively mono-energetic, and 
(ii) significant emission beyond this limit implies the emit- 
ting region is moving along the line of sight at relativistic 
speeds. Observations of the April 2011 flare suggest that, 
regarding (ii), a lower limit on the Doppler factor of 5 J> few 
is required. As for (i), the observed SED suggests a pile-up 
distribution that is not yet effectively mono-energetic. Inter- 



estingly, as Buehler et al. (20121 point out and we confirm, 



their data are not consistent with typical shock accelera- 



tion models, which usually produce p 2 (e.g. Kirk et al. 
2000 1. Instead, we suggest their observations are consistent 



with harder distributions (p < 1) found in many magnetic 



reconnection model^ 




Romanova & Lovelace 1992 


Zenitani 


& Hoshino|2001 2007 


. Notably, some reconnection models 



3 Note that |de Gouveia dal Pino fc Lazarian] l |2005| construct a 
reconnection model that produces p = 2 to 2.5 power-law indices, 



though O'C. Drury 1 2012 I argues against their analysis. 
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even predict a p = 1 electron distribution (Larrabee et al. 
20031 lO'C. Drury||2012l). Thus, the April 2011 flare SEDs 



may be consistent with emitting particles accelerated in a 
reconnection region that are undergoing bulk relativistic mo- 
tion with a Doppler factor of a few or more. 



5 A MINUET STATISTICAL MODEL 

To illustrate how reconnection minijets relate to the high 
energy nebula flux and variability, we construct a toy model 
that produces statistical predictions about the high energy 
nebula light curve. We postulate that reconnection minijets 
are random independent events in the nebula with an asso- 
ciated average reconnection event rate, n r , and are there- 
fore described by Poisson statistics. A significant simplify- 
ing assumption we make is in presuming that the probabil- 
ity density functions (PDFs) for the intrinsic reconnection 
emission region timescale, t', unbeamed intrinsic flux, /', 
and Lorentz factor, T, are narrow enough to be treated as 
Dirac delta probability densities. Thus, because r', /', and 
F are constants, the statistics of the random variables r (ob- 
served timescale) and / (observed flux), are determined by 
the minijet Doppler factor, 8, itself a random variable. An- 
other significant simplifying assumption we make is that the 
reconnection outflows have an isotropic angular distribution. 

The following discussion divides into two sections: |5.1| 
covers the statistics of individual minijets, and §5.2| develops 
time series statistics relevant to the nebular light curve as a 
whole. 



5.1 Individual minijet statistics 

Define a spherical coordinate system (r, <j>, 6) with the z-axis 
(9 — 0) pointing along the line of sight so that the viewing 
angle, 9, of any given jet is equal to the coordinate 6 associ- 
ated with its trajectory. Thus, the PDF p(5) is a function of 
the random variable 9 alone. Assuming the minijet angular 
distribution is isotropic so that p(5(9))dS = d(cos9), then 
from the definition of the Doppler factor, 



d ( cos^) = id(l-i) 



prs 2 



-dS. 



Therefore, the Doppler factor PDF is 

1 



p(8)dS 



I3YP 



dd. 



(13) 



(14) 



where <5 min = (r)" 1 and <5 max = (r-^r 2 - l))" 1 . We note 
that (5 m i n may effectively take on higher values, changing the 
normalization as well, so different values for cases where the 
debeamed flares gamma-ray flux is below a telescope's . 



Using equation ( 14 1 we now calculate observable quan- 



tities such as minijet timescale and flux distributions. Re- 
garding timescales, we find n(r)dr, which is the number of 
minijets that activate per unit time whose observed duration 
is between r and r + dr. This can be calculated by substi- 



tuting S = r'/r into equation (141 and including a factor of 



n(r)dr = . fur 



<t <Tt' 



(15) 



Hence, there is equal probability that short or long duration 




Figure 4. Cartoon schematic of reconnection sites in the nebula 
as viewed from above the toroidal plane (defined by the pulsar 
spin axis). As shown in the upper right corner box, each reconnec- 
tion site consists of plasma inflows into the reconnecting region 
and twin relativistic outflows ( "minijets" ) with some Lorentz fac- 
tor r of order few. Also shown in the schematic is a minijet di- 
rected toward the observer, which causes a flare since its emission 
is more highly beamed compared to its off-axis counterparts. 



minijets will be observed. To obtain the flux distribution, we 
assume the intrinsic gamma-ray flux of each mini-jet, /' is 
Doppler boosted such that / = 5 q f', where q — 3 + a (F e oc 
e~ a ) for moving components, and for bolometric flux, q — 4 
(Lind & Blandford||1985| |Jester||2008|). We now substitute 

to find the 



the Doppler boosting formula into equation 
minijet flux probability density 



P(f)df = 



qpTf 



df. 



(16) 



As expected, this formula is identical to that used in the cal- 
culation of luminosity functions for Doppler beamed sources 
(e.g. eqn. 2 of Urry & Shafcr 19841. This equation implies 
that p(f) oc / for 9 > 1, which means the minijet flux 
average, averaged over different flares, is dominated by rare 
bright flares. 



5.2 Minijet time series statistics 

In this subsection we derive quantities relevant to the entire 
high energy nebula light curve. To do this we analyze the 
statistics related to the random variable F na h, representing 
the total nebular high energy flux. First, we derive the sta- 
tistical behavior of the random variable k, or the number of 
flares active at any given moment, and then apply this to 
the analysis of F nc h- Finally, assuming all nebular variability 
is due to minijets, we derive the nebular light curve's power 
spectrum. 

The statistical behavior of the light curve depends sen- 
sitively on A, the average number of reconnection events that 
are active at any given time. To derive A, we first note that 
the differential flare rate, dn/dS, immediately follows from 
equation (141 as dn/dS = n r p(5) (recall n r = average nebu- 
lar reconnection rate). Thus, the differential overlap number 
is dX = rdn (recall r is the observed flare duration), since 



rdn is the differential number of flares with Doppler fac- 
tor 5 that are activate within one observed flare duration r. 
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To calculate the average overlap number over an interval of Where of t = ((Fr) — {Fr) 



2\l/2 



Doppler factors from 8 n 
A : 



to S n 



we integrate to obtain 



is the standard deviation 

2\ i o / rp \3\ 



r'(«5 



2 

max 



5 2 - ) 

-'mm/ 



n r r 



(for 8 

max min J 



(17) 



Since flaring consists of a random process of unrelated 
events governed by an average overlap number A within the 
nebula, the probability that k flares are active/overlapping 
at any given moment is governed by the Poisson distribution, 
P(k) = \ k e~ x /k\. We have verified the applicability of the 



Poisson distribution and our calculation of A (eqn. 171 to 
our model via Monte Carlo simulations with a variety of 
parameters values. 

The duty cycle of minijets (i.e. the fraction of time when 
the nebula contains one or more active minijet) can now be 
calculated. The duty cycle of flares with <5 m i n < 5 < (5 max is 
simply the probability that at a time t, one or more flares 
(k 1) are active: /duty = P(k ^ 1) = 1 — e~ A . In the non- 
overlapping limit, A< 1, the duty cycle reduces to /duty ~ 
A. 

We may now analyze the statistics of F ncCl , the random 
variable representing the total high energy flux of the neb- 
ula. To this end, we introduce the random variable for the 
total minijet flux due to reconnection outflows, F r , and the 
constant nebula flux, -Fconst, possibly originating from the 
pulsar wind termination shock; the random variables F ne h 
and F r are related by 



Fr — F nt 



(18) 



Note that if only one minijet were active at time to, then 
F r {to) — f. In order to derive the statistical behavior of 
the directly observable -F ne b, we analyze the PDF of F r , 
p(F r ), where p(F r )dF r is the probability that, at any given 
time, the high energy flux due to reconnection minijets is 
between F r and F r +dF r . Each minijet is presumed to have a 
square pulse profile. Importantly, our analysis assumes that 
p(F r ) and its parameters remain constant over time intervals 
well in excess of the light curve's autocorrelation time, T a u to . 
If this assumption holds, then the moments of p(F r ) may 
be compared with the observed moments of the light curve 
extracted from data spanning time intervals larger than the 
observed r au to. 

The exact moment generating function of p(F r ) can be 
found, which contains the same information as p(F r ) (see 
appendix for its derivation). While any moment of p(F r ) 
may be calculated, we provide the first three moments here 
for /' = 1: 



/ jp \ _ n r T ( fl-2/q _ fl-2/q 
\- rr / — („ q\ «-p W max -'min 



(«-2)j9T 



9-2 



m 1/2 (f^ /q - fi- n 2/q 



1/2 



<F r ) (29-2)1/=" inrT ,y,2 

(2g - 2 )»/» {Prr 2 {f^ q - f^ ,q ) 

3 ? -2 , „_Al/2/ , 2-2/ 9 _ f 2-2/ q \ 3 / 2 
(I It I ) Jmax ./min 



, (19) 



and 7i = cr" 3 ((F 3 ) - 3 (F r ) (F 2 ) + 2 <F r ) 3 ) is the skew- 
ness; the brackets () represent a time average. As a check 
against our analytical work here, we have spot checked equa- 



tions ( 19 1 with Monte Carlo simulations using a variety of 



parameters and found them to be consistent with one an- 
other. Note that comparing the above moments to an ob- 
served light curve's moments allows for one to test whether 
the light curve is generated by beamed minijets per our 
model. To better discuss the significance of equations uXsh , 
we assume / ma x > /min, T > 1, and / max w (2P) 9 , so that 



equations ( 19 1 can be approximated as 
2 q - 2 n r r'V>-' A 

{Fr) 

7i 



9 - 

(2g- 



2 

4)pV 2 



{{2q-2)n r T'y/ 2 
2(2g-2) 3/2 P 3/2 



(20) 



(21) 



(22) 



(3g-2) (rvr')V2' 

The variation about the average as expressed by the rela- 
tive standard deviation {of t I (F r )) can easily take on small 
or large values depending on how n r r' compares with T 3 . 
Because 71 depends on V, n r , and r' in the same way as 
a I (F r ) in this approximation, their ratio depends only on 
the beaming index and is expected to be of order unity: 

o*, „ (2g~4)(3g-2) 
(F r ) 7l ~ 2(2g - 2)2 1 ' 

where this ratio for reasonable values of the beaming factor 
(q > 3) ranges from cr F J((F r ) 71) ~ 0.44 to 0.75. 

Another important statistical representation of the light 
curve is the power spectrum, P(v), which measures its vari- 
ability on different timescales. We compute P(y) for light 
curves composed of square pulses, exponential pulses (i.e. 
zero rise time and exponential decay), and Gaussian pulses, 
all with a timescale of r = t'/5. Since individual pulses are 
generated by a stochastic process, the phase information of 
the Fourier transform is unimportant, implying the power 
spectra for an individual pulse, Pi{v), can be used to gen- 
erate the light curve power spectrum thusly: 



P(u) = 



n{T)Pi{v)dr, 



(24) 



where n(r) is found in equation (151. For the square pulse 



and the exponential pulse, the single pulse power spectra 
are Pp oc sin (ttvt) 2 /v 2 and P 7 cxp oc (1 + 4tt 2 t 2 // 2 )- 1 re- 
spectively (do to excessive length we do not include the 
spectrum for Gaussian pulses). For the total power spec- 
trum, when v <C T"mix ~ (2rr') _1 , the power spectrum is 
constant (i.e. white noise) since the frequency dependence 
drops out of the integrand. For v ^> r mm m 2F/t' , only the 
power spectrum for square and exponential pulses goes as 
P(v) oc v~ 2 (for the square pulses, there is oscillatory be- 
havior on top of the v~ 2 decay). This can be summarized as 
follows: 

{const. for v < (2rr') _1 

transition to u~ 2 for (2rr') _1 < v < 2T/t' 
v~ 2 for v > 2r/r' 

(25) 

Unlike the square and exponential pulses, Gaussian pulses 



© 0000 RAS, MNRAS 000, [Tpj] 



Crab flares as minijets 9 



Exponential and Gaussian pulses 




10"' 10" 
Square Pulses 



P(v) « v + oscillations 




10"' 10" 
v (days -1 , for x min = 4 days and r=3) 

Figure 5. Power spectra for exponential, Gaussian, and square 
pulses. The exponential and Gaussian power spectra (top) dis- 
play a marked difference, in that the exponential power spectrum 
decays as a power-lay (oc v~ 2 ) and the Gaussian power spectrum 
falls off more rapidly. The square pulse power spectrum (bottom) 
is similar to the exponential one since it goes as u~ 2 , albeit with 
superimposed sinusoidal oscillations. The gray data represents the 
square of the discrete Fourier transform (DFT) of the Monte Carlo 
generated light curve via square pulses. The normalization of the 
square pulse P(v) curve in the bottom panel is adjusted by hand 
for easy comparison with the simulated DFT (the exponential 
and Gaussian P(v) curves use the same normalization). 



do not display power-law, they drop off faster (see fig [5| . 
While we have found analytical expressions for P(v) for 
some beaming factors (e.g. q = 4), we do not reproduce 
them here due to their excessive length. However, the so- 
lutions to equation ( |24| ) are plotted in figure [5] As a check 
against our analytical work, we include the discrete Fourier 
transform of a Monte Carlo simulated light curve composed 
of square pulses, which clearly follows the corresponding cal- 
culated power spectrum. 



5.3 Application of Statistical Model to Crab 
Nebula 

Our model may be divided into two categories: the non- 
overlapping regime and the overlapping regime. In the non- 
overlapping regime, minijets do not temporally overlap and 
a significant constant emission component dominates, i.e. 
A « 1 and F const / (F r ) 3> 1. The opposite is true for the 
overlapping regime, where A 3> 1 and .Fconst/ (F r ) < 1. 
In the non-overlapping regime, each observed flare is due 
to an individual minijet's flux, /, hence individual mini- 
jet statistics are directly observable. In this approximation 
fconst/ (F r ) 2> 1, so that Fconst can be related to an observ- 
able quantity via (F neD ) ~ Fconst- Therefore, we can write 
/ « F flarc - (Fncb), where F flaro is the nebula flux (F ne b) 



> 

2 
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Figure 6. Ten year simulated Crab nebula light curve. The 
"April 2011 type flares" represent flares with increases of ~ 
30 over the nebular average as found in |Buehler et al.| | |2012| >. 
"September type 2010" flares represent flux increases by ~ 5 sim- 
ilar to the Crab September 2010 flare jAbdo et al.|2011| |Tavani| 
|et al.||201l) . The timescale assumes the shortest flare durations 
are 4 days, thus r' ~ 2F X 4 days, so if T = 3, then r' = 24 days. 
The simulated light curve was binned into 4 day intervals. For 
easy comparison to data, we add a steady component of > 100 
MeV flux Fconst = 6 X 10 — 7 s — 'cm -2 , consistent with the Crab 
Nebula's average flux as measured by Fermi/LAT | |Abdo et al.| 
|2011| . The model parameters are (5 m im <5 m ax, Q> n n T ' > F, /') 



(0.56,5.83,4,0.1875 days -1 , 24 days, 3, 0.15 X 10~ 7 s 
with a corresponding flare overlap number of A = 1.25. 



during an isolated flare. If we insert this expression into 

(26) 



equation ( 16 1 we obtain 



p(F na rc)dFfl a r c 

oc (F n a 



9 + 1 

(Fneb)) * dF R 

arc 



again, where F narc = F I1CD when a single flare is active. 

A separate observable in this regime is the minijet 



timescale distribution (eqn. 15 1 which states that flares will 
be equally distributed across the allowed timescales. 

For the overlapping regime, the above analysis is 
complicated by the difficulty of associating a single flare 
with a single minijet. However, applicable to both regimes 
is the light curve power spectrum (eqn. 251, which for 
some pulse profiles is oc v~ 2 in the short variability region. 



However, this has been measured by Buehler et al. (20121 
as being oc i>~ 1 . Such a measurement could be explained if 
the transition region where P v evolves from being constant 



to being oc 



is dominating that measurement. Other 



19 1, such 



observables are the light curves moments (eqns. 
as the ratio of variability to skewness (eqn. 231. However, 
moments such as the relative standard deviation, which 
measures the light curve variability, cannot be directly 
compared with observational data since observed light 
curves smooth over short timescale variability because of 
time binning and the breaks in observational coverage. 
As an example light curve, figure [6] represents of the 
non-overlapping regime of our model that qualitatively 
reproduces the Crab light curve. The parameters used 
to generate figure [6] are (<5 m i n , <W>c, q, n r , r', T, /') = 
(0.56,5.83,4,0.1875 days" 1 , 24 days, 3, 0.15 x 
10 -7 s -1 cm -2 ), which imply a flare overlap number of 
A = 1.25. 

For either regime, the energy in PWNe reconnection 
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outflows must not exceed the total energy budget set by 
the spin-down power of the pulsar. For example, if the 
non-overlapping model is true and approximately one flare 
per year is observed with a viewing ang lc of < r~\ then 
n r ~ 4F 2 yr -1 . For the September 2010 flare parameters, 
this implies the minijet power expended is 

P r = n r E r = 4.0 x lO^Lae.eis.srv/r 1 erg s"\ (27) 

which is much less than the Crab pulsar's spin-down power, 
Pspin « 5 x 10 38 erg s _1 . 

In applying our statistical minijet model to the Crab 
Nebula, its simplifying assumptions that (a) the minijets' 
intrinsic parameters are the same and that (b) the minijet 
directions are isotropically distributed are both open to crit- 
icism. Assumption (a) is challenged by the most luminous 
flare being longer in duration (~ 9 days) than the Septem- 
ber 2010 flare (~ 4 days), since both assumption (a) and the 
Doppler transformations indicate the April 2011 flare should 
be shorter. However, more flare observations are necessary 
before any firm conclusions are made regarding a correla- 
tion (or lack thereof) between observed flare luminosity and 
duration. Regarding assumption (b), the clear toroidal mor- 
phology of the Crab nebula (e.g. Weisskopf et al.|2000 1 , con- 
sistent with pulsar models with a toroidal outflow containing 
a toroidally dominated magnetic field and a large-scale cur- 
rent sheet, combine to make the assumption of an isotropic 
angular distribution of minijets questionable. In the PWN 



split monopole model or the striped wind one ( Coroniti 1990 
Bogovalov| 19991 ), reconnection in the implied current sheets 



would only produce minijets in the plane of the torus, ren- 
dering them unobservable since the pulsar spin axis is at an 
angle of ~ 60° to our line of sight ( Weisskopf et al. 2000 1 . For 



this reason, we have assumed the minijets are produced in 
an turbulent isotropic region of the nebula. Nonetheless, the 
large-scale anisotropic morphology and structure of the mag- 
netic field suggests that future work on this model should 
take into account at least some degree of anisotropy. 



such distributions can easily form a pile-up near the radi- 
ation reaction limit implied by MHD, effectively becoming 



mono-energetic. The April 2011 flare SEDs (Buehler et al 



2012 1 are inconsistent with shock acceleration and are in- 



stead consistent with a hard electron distribution that may 
contain a not yet effectively mono-energetic pile-up. Fur- 
thermore, because the observed location of the SED peak 
is greater than the peak predicted by the synchrotron radi- 
ation reaction limit, the April 2011 flare's emitting region 
Doppler factor is > few. 

The predictions of our statistical minijet model can be 
summarized as follows: 

• When minijets do not temporally overlap one another, 
the probability density function for the nebula's high energy 
flux during a flare, Ffl are , is 



P(Fflarc) ~ COnSt (ffl a 



CFneb})" 



(28) 



• The first three moments of the light curve may be com- 
pared with our theoretically calculated moments in equa- 



tions (191, and any higher degree theoretical moments may 



be easily calculated using the method described in the ap- 
pendix. 



• The light curve power spectrum (eqn. 25 1 is constant 
("white noise") for v <C (Ft') -1 , and goes as P(y) oc v~ 2 
for v » T/r. 

Unlike the standard model for Crab Nebula non- 
thermal emission, wherein particles are accelerated by the 



pulsar wind termination shock (Kennel & Coroniti 1984), 



we have suggested here that magnetic reconnection may 
play an important or even dominant role in particle accel- 
eration. Further research will consider whether this recon- 
nection model can explain both the steady nebula emission 
and the flaring and therefore preclude the need for shock 
emission altogether. Our statistical model may also apply 
for AGNs that exhibit several minute TeV variability. 



6 CONCLUSIONS 

We have constructed a statistical model of Crab Nebula high 
energy variability by assuming that magnetic reconnection 
sites throughout the nebula are activated randomly, and 
once activated, display emission characteristics controlled 
by the reconnection outflow Doppler factor. At each site, a 
magnetically dominated reconnection region launches twin 
relativistic outflows along a randomly aligned axis. GeV 
flares are observed when, by chance, a relativistic outflow is 
aligned with the line of sight and is thus Doppler-boosted. 
The observed flares' unusually short durations and high lu- 
minosities suggest the emitting plasma is indeed moving to- 
ward Earth at relativistic speeds. 

The flares' SEDs contain information about the particle 
acceleration process that suggests non-thermal particles are 
generated in reconnection regions, rather than by shocks, 
and are undergoing bulk relativistic motion along the line of 
sight. Models of reconnection particle acceleration tend to 
create a hard power-law with an index close to p ~ 1 ( |Ro-| 
manova fc Lovelace|[T992l |Zenitani fc Hoshmol[200l| |2007[ 



Larrabee et al.|2003||O'C. Drury|2012[ ). We have shown that 



ACKNOWLEDGMENTS 

We thank the Fermi/LAT team for making data concerning 
the April 2011 flare public and Konstantinos Gourgouliatos 
for many helpful discussions. 



REFERENCES 

Abdo et al. 2010, ApJ, 708, 1254 
Abdo et al. 2011, Science, 331, 739 

Balbo M., Walter R., Ferrigno C, Bordas P., 2011, A&A, 
527, L4+ 

Bednarek W., Idee W., 2011, MNRAS, 414, 2229 
Begelman M. C, 1998, ApJ, 493, 291 

Blackman E. G., Field G. B., 1994, Physical Review Let- 
ters, 72, 494 
Bogovalov S. V., 1999, A&A, 349, 1017 
Buehler et al. 2012, ApJ, 749, 26 

Bulmer M. G., 1979, Principles of Statistics. New York, 

Dover Publications. 
Cerutti B., Uzdensky D. A., Begelman M. C, 2012, ApJ, 

746, 148 

Coroniti F. V., 1990, ApJ, 349, 538 



© 0000 RAS, MNRAS 000, [Tpj] 



Crab flares as minijets 11 



de Gouveia dal Pino E. M., Lazarian A., 2005, A&A, 441, 
845 

de Jager O. O, Harding A. K., Michelson P. F., Nel H. I., 
Nolan P. L., Sreekumar P., Thompson D. J., 1996, ApJ, 
457, 253 

Espinoza C. M., Jordan C, Stappers B. W., Lyne A. G., 
Weltevrede P., Cognard I., Theureau G., 2010, The As- 
tronomer's Telegram, 2889, 1 
Gallant Y. A., Arons J., 1994, ApJ, 435, 230 
Giannios D., Uzdensky D. A., Begelman M. C, 2009, MN- 
RAS, 395, L29 

Giannios D., Uzdensky D. A., Begelman M. C, 2010, MN- 
RAS, pp 30-+ 

Hester J. J., Mori K., Burrows D., Gallagher J. S., Graham 

J. R., Halverson M., Kader A., Michel F. C, Scowen P., 

2002, ApJ, 577, L49 
Hester et al. 1995, ApJ, 448, 240 
Jester S., 2008, MNRAS, 389, 1507 
Kennel C. F., Coroniti F. V., 1984, ApJ, 283, 694 
Kirk J. G., Guthmann A. W., Gallant Y. A., Achterberg 

A., 2000, ApJ, 542, 235 
Kirk J. G., Rieger F. M., Mastichiadis A., 1998, A&A, 333, 

452 

Kirsch et al. 2005, in O. H. W. Siegmund ed., Society of 
Photo-Optical Instrumentation Engineers (SPIE) Confer- 
ence Series Vol. 5898 of Society of Photo-Optical Instru- 
mentation Engineers (SPIE) Conference Series, Crab: the 
standard x-ray candle with all (modern) x-ray satellites, 
pp 22-33 

Komissarov S. S., Lyubarsky Y. E., 2004, MNRAS, 349, 
779 

Komissarov S. S., Lyutikov M., 2011, MNRAS, 414, 2017 
Larrabee D. A., Lovelace R. V. E., Romanova M. M., 2003, 

ApJ, 586, 72 
Lind K. R., Blandford R. D., 1985, ApJ, 295, 358 
Lyubarsky Y., Kirk J. C, 2001, ApJ, 547, 437 
Lyubarsky Y. E., 2005, MNRAS, 358, 113 
Lyutikov M., 2006, MNRAS, 369, L5 
Lyutikov M., 2010, MNRAS, 405, 1809 
Lyutikov M., Balsara D., Matthews C, 2011, ArXiv e- 

prints: 1109.1204 
Lyutikov M., Blackman E. C, 2001, MNRAS, 321, 177 
Lyutikov M., Blandford R., 2003, ArXiv Astrophysics e- 

prints 

Lyutikov M., Uzdensky D., 2003, ApJ, 589, 893 
McKinney J. C, Uzdensky D. A., 2010, ArXiv e-prints: 
1011.1904 

Melrose D. B., 1980, Plasma astrophysics. Nonthermal pro- 
cesses in diffuse magnetized plasmas - Vol.1: The emission, 
absorption and transfer of waves in plasmas; Vol.2: Astro- 
physical applications 
Michel F. C, 1973, ApJ, 180, L133+ 
O'C. Drury L., 2012, ArXiv e-prints: 1201.6612 
Priest E., Forbes T., 2000, Magnetic Reconnection 
Rees M. J., Gunn J. E., 1974, MNRAS, 167, 1 
Rees M. J., Meszaros P., 1994, ApJ, 430, L93 
Romanova M. M., Lovelace R. V. E., 1992, A&A, 262, 26 
Schlickeiser R., 1984, A&A, 136, 227 

Spada M., Ghisellini G., Lazzati D., Celotti A., 2001, MN- 
RAS, 325, 1559 
Striani et al. 2011, ApJ, 741, L5+ 
Tavani et al. 2011, Science, 331, 736 

© 0000 RAS, MNRAS 000,|ilfl3l 



Trimble V., 1983, Reviews of Modern Physics, 55, 511 
Urry C. M., Shafer R. A., 1984, ApJ, 280, 569 
Uzdensky D. A., 2011, Space Sci. Rev., pp 101 — h 
Uzdensky D. A., Cerutti B., Begelman M. C, 2011, ApJ, 
737, L40 

Weisskopf et al. 2000, ApJ, 536, L81 

Wilson-Hodge et al. 2011, ApJ, 727, L40+ 

Yuan Q., Yin P.-F., Wu X.-F., Bi X.-J., Liu S., Zhang B., 

2011, ApJ, 730, L15+ 
Zenitani S., Hoshino M., 2001, ApJ, 562, L63 
Zenitani S., Hoshino M., 2007, ApJ, 670, 702 



12 E. Clausen-Brown and M. Lyutikov 



APPENDIX A: CALCULATING THE MOMENTS OF p(F R ) 

Here we exploit a special property of moment generating functions (MGFs, e.g. Bulmer|1979 1 that allows us to calculate the 
MGF for the probability density function (PDF), p(F r ). First we state the expression for the probability that the observed 
minijet flux will be between F r and F r + dF r at any given time: 



A e 



p(F r )dF r =^2P(k)p(F k )dF r = S(F r )e- x dF r + ^ — — P (F k )dF, 



(Al) 



where for compactness we have replaced F r \k, or "F r given k", with F k . It is important to note that p(F r ) is the PDF for the 
total summed minijet flux at any given timeslice, not the PDF that governs the flux of indiviudal flares, which is described 
by equation ( 16 1. 

The MGF for a random variable Y that follows the PDF p(Y), is defined as My ~ exp (Yx)p(Y)dY , which can be 
Taylor expanded to produce 



M Y = 1 + (Y) x + - (Y 2 ) . 



Y J )x 



(A2) 



where i is a dummy variable such that any moment of p(Y) may be calculated via (Y 1 ) — d 3 My /dx J \ x =o- Applying the 
definition of My to p(F r ) we find 



M Fr {x) = (e^) = e- x 



fc! 



(A3) 



where Mf,. is the MGF associated with conditional PDF, p(F k ). 

The special property of MGF s we use here is that for a random variable Ft that is the sum of k random variables, 
F k = ELi Aj the MGF is simply (|Bulmer||l979 \ 



Si' 



(A4) 



where M/ 4 is the moment generating function for the random variable fi, the i* minijet flux of the k minijets active at time 
t (the more debeamed minijet in each reconnection outflow pair is not taken into account as it contributes little to the overall 
flux). Because the PDFs for all minijet random variables, fi, are the same, then 



M F , = M 



(A5) 



Mf is easily obtained from calculating the moments of p(Fk=i), which is equation (16 1 with an extra factor of r and a 
correspondingly different normalization. This extra factor of r is necessary because longer flares are more likely to be active 
at a time t compared with shorter duration flares, while equation ( 16 1 does not take a flare's time dependence into account. 
We now calculate (ff): 



FY = 



F£ =1 p(F k=1 )dF k=1 = A 



/min 



fp(f)rdf 



/min 



tj f 2 /i _ fit fi 

JmaxJ m i n JmaxJ m [ n 
f2/q _ f2/q 
J max J min 



2/3 fj 



(A6) 



where / is the single flare flux and A is a normalization constant. Now Mp x can be obtained using equation \A2\ , such that 
AIf 1 — E which can be inserted into equation (A5l to produce 



/ oo 

M Fk =M Fl = [ZjM 
\j=o J ' 



Hence, the total MGF associated with p(F r ) is 



= e 



fc=l \j=0 J J 



(A7) 



(A8) 



To find the first three moments of Mf t , we calculate the first three moments of Mf^ by expanding Mf h to third order in the 
dummy variable x: 



Mf, = M 



l + k(F 1 )x+[~ (F, 2 ) + 



P 3\ , fc(fc-l) 



.. , fc(fc-l)(fc-2) 



(A9) 
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Therefore, we can now read off the first three moments of the p(Ffc) distribution 

{F k ) = k{F x ) 
(-Ffe ) = k (F 2 ) + k(k - 
(F^) = fc(F 3 )+3fc(fc 

Now, note that 



1) (Fi) 

■l){F 1 )(Ff) + k{k-l){k-2){F 1 ) 



Fi ) = e 



E 



k\ 



(A10) 



(All) 



Evaluating the infinite sum in equation (All \ for k = 1, 2, and 3 leads us to the first three moments of the p(F r ) distribution, 



(F r ) = \(F 1 ) 
(F r 2 ) = \(F?)+\ 2 (Fi) 2 
<F r 3 ) = A (Ff) + 3A 2 (Fx) (F 2 ) + A 3 (Fi) 3 



(A12) 



Calculating higher order moments is stra ight forward. In general, the j th moment of p(F r ) can be calculated by Taylor 
expanding the expressions for MF k (eqn. A9l, reading off the j th order term in said expansion, and using this term in 
evaluating the infinite sum for (F;?) (eqn.[All|. 
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